home *** CD-ROM | disk | FTP | other *** search
/ AGA Toolkit '97 / The AGA Toolkit '97.iso / miscellaneous / science / maths / calc / source / commath.c < prev    next >
Encoding:
C/C++ Source or Header  |  1996-09-07  |  9.4 KB  |  604 lines
  1. /*
  2.  * Copyright (c) 1993 David I. Bell
  3.  * Permission is granted to use, distribute, or modify this source,
  4.  * provided that this copyright notice remains intact.
  5.  *
  6.  * Extended precision complex arithmetic primitive routines
  7.  */
  8.  
  9. #include "cmath.h"
  10.  
  11.  
  12. COMPLEX _czero_ =        { &_qzero_, &_qzero_, 1 };
  13. COMPLEX _cone_ =        { &_qone_, &_qzero_, 1 };
  14. COMPLEX _conei_ =        { &_qzero_, &_qone_, 1 };
  15.  
  16. static COMPLEX _cnegone_ =    { &_qnegone_, &_qzero_, 1 };
  17.  
  18.  
  19. /*
  20.  * Free list for complex numbers.
  21.  */
  22. static FREELIST freelist = {
  23.     sizeof(COMPLEX),    /* size of an item */
  24.     100            /* number of free items to keep */
  25. };
  26.  
  27.  
  28. /*
  29.  * Add two complex numbers.
  30.  */
  31. COMPLEX *
  32. cadd(c1, c2)
  33.     COMPLEX *c1, *c2;
  34. {
  35.     COMPLEX *r;
  36.  
  37.     if (ciszero(c1))
  38.         return clink(c2);
  39.     if (ciszero(c2))
  40.         return clink(c1);
  41.     r = comalloc();
  42.     if (!qiszero(c1->real) || !qiszero(c2->real))
  43.         r->real = qadd(c1->real, c2->real);
  44.     if (!qiszero(c1->imag) || !qiszero(c2->imag))
  45.         r->imag = qadd(c1->imag, c2->imag);
  46.     return r;
  47. }
  48.  
  49.  
  50. /*
  51.  * Subtract two complex numbers.
  52.  */
  53. COMPLEX *
  54. csub(c1, c2)
  55.     COMPLEX *c1, *c2;
  56. {
  57.     COMPLEX *r;
  58.  
  59.     if ((c1->real == c2->real) && (c1->imag == c2->imag))
  60.         return clink(&_czero_);
  61.     if (ciszero(c2))
  62.         return clink(c1);
  63.     r = comalloc();
  64.     if (!qiszero(c1->real) || !qiszero(c2->real))
  65.         r->real = qsub(c1->real, c2->real);
  66.     if (!qiszero(c1->imag) || !qiszero(c2->imag))
  67.         r->imag = qsub(c1->imag, c2->imag);
  68.     return r;
  69. }
  70.  
  71.  
  72. /*
  73.  * Multiply two complex numbers.
  74.  * This saves one multiplication over the obvious algorithm by
  75.  * trading it for several extra additions, as follows.  Let
  76.  *    q1 = (a + b) * (c + d)
  77.  *    q2 = a * c
  78.  *    q3 = b * d
  79.  * Then (a+bi) * (c+di) = (q2 - q3) + (q1 - q2 - q3)i.
  80.  */
  81. COMPLEX *
  82. cmul(c1, c2)
  83.     COMPLEX *c1, *c2;
  84. {
  85.     COMPLEX *r;
  86.     NUMBER *q1, *q2, *q3, *q4;
  87.  
  88.     if (ciszero(c1) || ciszero(c2))
  89.         return clink(&_czero_);
  90.     if (cisone(c1))
  91.         return clink(c2);
  92.     if (cisone(c2))
  93.         return clink(c1);
  94.     if (cisreal(c2))
  95.         return cmulq(c1, c2->real);
  96.     if (cisreal(c1))
  97.         return cmulq(c2, c1->real);
  98.     /*
  99.      * Need to do the full calculation.
  100.      */
  101.     r = comalloc();
  102.     q2 = qadd(c1->real, c1->imag);
  103.     q3 = qadd(c2->real, c2->imag);
  104.     q1 = qmul(q2, q3);
  105.     qfree(q2);
  106.     qfree(q3);
  107.     q2 = qmul(c1->real, c2->real);
  108.     q3 = qmul(c1->imag, c2->imag);
  109.     q4 = qadd(q2, q3);
  110.     r->real = qsub(q2, q3);
  111.     r->imag = qsub(q1, q4);
  112.     qfree(q1);
  113.     qfree(q2);
  114.     qfree(q3);
  115.     qfree(q4);
  116.     return r;
  117. }
  118.  
  119.  
  120. /*
  121.  * Square a complex number.
  122.  */
  123. COMPLEX *
  124. csquare(c)
  125.     COMPLEX *c;
  126. {
  127.     COMPLEX *r;
  128.     NUMBER *q1, *q2;
  129.  
  130.     if (ciszero(c))
  131.         return clink(&_czero_);
  132.     if (cisrunit(c))
  133.         return clink(&_cone_);
  134.     if (cisiunit(c))
  135.         return clink(&_cnegone_);
  136.     r = comalloc();
  137.     if (cisreal(c)) {
  138.         r->real = qsquare(c->real);
  139.         return r;
  140.     }
  141.     if (cisimag(c)) {
  142.         q1 = qsquare(c->imag);
  143.         r->real = qneg(q1);
  144.         qfree(q1);
  145.         return r;
  146.     }
  147.     q1 = qsquare(c->real);
  148.     q2 = qsquare(c->imag);
  149.     r->real = qsub(q1, q2);
  150.     qfree(q1);
  151.     qfree(q2);
  152.     q1 = qmul(c->real, c->imag);
  153.     r->imag = qscale(q1, 1L);
  154.     qfree(q1);
  155.     return r;
  156. }
  157.  
  158.  
  159. /*
  160.  * Divide two complex numbers.
  161.  */
  162. COMPLEX *
  163. cdiv(c1, c2)
  164.     COMPLEX *c1, *c2;
  165. {
  166.     COMPLEX *r;
  167.     NUMBER *q1, *q2, *q3, *den;
  168.  
  169.     if (ciszero(c2))
  170.         math_error("Division by zero");
  171.     if ((c1->real == c2->real) && (c1->imag == c2->imag))
  172.         return clink(&_cone_);
  173.     r = comalloc();
  174.     if (cisreal(c1) && cisreal(c2)) {
  175.         r->real = qdiv(c1->real, c2->real);
  176.         return r;
  177.     }
  178.     if (cisimag(c1) && cisimag(c2)) {
  179.         r->real = qdiv(c1->imag, c2->imag);
  180.         return r;
  181.     }
  182.     if (cisimag(c1) && cisreal(c2)) {
  183.         r->imag = qdiv(c1->imag, c2->real);
  184.         return r;
  185.     }
  186.     if (cisreal(c1) && cisimag(c2)) {
  187.         q1 = qdiv(c1->real, c2->imag);
  188.         r->imag = qneg(q1);
  189.         qfree(q1);
  190.         return r;
  191.     }
  192.     if (cisreal(c2)) {
  193.         r->real = qdiv(c1->real, c2->real);
  194.         r->imag = qdiv(c1->imag, c2->real);
  195.         return r;
  196.     }
  197.     q1 = qsquare(c2->real);
  198.     q2 = qsquare(c2->imag);
  199.     den = qadd(q1, q2);
  200.     qfree(q1);
  201.     qfree(q2);
  202.     q1 = qmul(c1->real, c2->real);
  203.     q2 = qmul(c1->imag, c2->imag);
  204.     q3 = qadd(q1, q2);
  205.     qfree(q1);
  206.     qfree(q2);
  207.     r->real = qdiv(q3, den);
  208.     qfree(q3);
  209.     q1 = qmul(c1->real, c2->imag);
  210.     q2 = qmul(c1->imag, c2->real);
  211.     q3 = qsub(q2, q1);
  212.     qfree(q1);
  213.     qfree(q2);
  214.     r->imag = qdiv(q3, den);
  215.     qfree(q3);
  216.     qfree(den);
  217.     return r;
  218. }
  219.  
  220.  
  221. /*
  222.  * Invert a complex number.
  223.  */
  224. COMPLEX *
  225. cinv(c)
  226.     COMPLEX *c;
  227. {
  228.     COMPLEX *r;
  229.     NUMBER *q1, *q2, *den;
  230.  
  231.     if (ciszero(c))
  232.         math_error("Inverting zero");
  233.     r = comalloc();
  234.     if (cisreal(c)) {
  235.         r->real = qinv(c->real);
  236.         return r;
  237.     }
  238.     if (cisimag(c)) {
  239.         q1 = qinv(c->imag);
  240.         r->imag = qneg(q1);
  241.         qfree(q1);
  242.         return r;
  243.     }
  244.     q1 = qsquare(c->real);
  245.     q2 = qsquare(c->imag);
  246.     den = qadd(q1, q2);
  247.     qfree(q1);
  248.     qfree(q2);
  249.     r->real = qdiv(c->real, den);
  250.     q1 = qdiv(c->imag, den);
  251.     r->imag = qneg(q1);
  252.     qfree(q1);
  253.     qfree(den);
  254.     return r;
  255. }
  256.  
  257.  
  258. /*
  259.  * Negate a complex number.
  260.  */
  261. COMPLEX *
  262. cneg(c)
  263.     COMPLEX *c;
  264. {
  265.     COMPLEX *r;
  266.  
  267.     if (ciszero(c))
  268.         return clink(&_czero_);
  269.     r = comalloc();
  270.     if (!qiszero(c->real))
  271.         r->real = qneg(c->real);
  272.     if (!qiszero(c->imag))
  273.         r->imag = qneg(c->imag);
  274.     return r;
  275. }
  276.  
  277.  
  278. /*
  279.  * Take the integer part of a complex number.
  280.  * This means take the integer part of both components.
  281.  */
  282. COMPLEX *
  283. cint(c)
  284.     COMPLEX *c;
  285. {
  286.     COMPLEX *r;
  287.  
  288.     if (cisint(c))
  289.         return clink(c);
  290.     r = comalloc();
  291.     r->real = qint(c->real);
  292.     r->imag = qint(c->imag);
  293.     return r;
  294. }
  295.  
  296.  
  297. /*
  298.  * Take the fractional part of a complex number.
  299.  * This means take the fractional part of both components.
  300.  */
  301. COMPLEX *
  302. cfrac(c)
  303.     COMPLEX *c;
  304. {
  305.     COMPLEX *r;
  306.  
  307.     if (cisint(c))
  308.         return clink(&_czero_);
  309.     r = comalloc();
  310.     r->real = qfrac(c->real);
  311.     r->imag = qfrac(c->imag);
  312.     return r;
  313. }
  314.  
  315.  
  316. /*
  317.  * Take the conjugate of a complex number.
  318.  * This negates the complex part.
  319.  */
  320. COMPLEX *
  321. cconj(c)
  322.     COMPLEX *c;
  323. {
  324.     COMPLEX *r;
  325.  
  326.     if (cisreal(c))
  327.         return clink(c);
  328.     r = comalloc();
  329.     if (!qiszero(c->real))
  330.         r->real = qlink(c->real);
  331.     r->imag = qneg(c->imag);
  332.     return r;
  333. }
  334.  
  335.  
  336. /*
  337.  * Return the real part of a complex number.
  338.  */
  339. COMPLEX *
  340. creal(c)
  341.     COMPLEX *c;
  342. {
  343.     COMPLEX *r;
  344.  
  345.     if (cisreal(c))
  346.         return clink(c);
  347.     r = comalloc();
  348.     if (!qiszero(c->real))
  349.         r->real = qlink(c->real);
  350.     return r;
  351. }
  352.  
  353.  
  354. /*
  355.  * Return the imaginary part of a complex number as a real.
  356.  */
  357. COMPLEX *
  358. cimag(c)
  359.     COMPLEX *c;
  360. {
  361.     COMPLEX *r;
  362.  
  363.     if (cisreal(c))
  364.         return clink(&_czero_);
  365.     r = comalloc();
  366.     r->real = qlink(c->imag);
  367.     return r;
  368. }
  369.  
  370.  
  371. /*
  372.  * Add a real number to a complex number.
  373.  */
  374. COMPLEX *
  375. caddq(c, q)
  376.     COMPLEX *c;
  377.     NUMBER *q;
  378. {
  379.     COMPLEX *r;
  380.  
  381.     if (qiszero(q))
  382.         return clink(c);
  383.     r = comalloc();
  384.     r->real = qadd(c->real, q);
  385.     r->imag = qlink(c->imag);
  386.     return r;
  387. }
  388.  
  389.  
  390. /*
  391.  * Subtract a real number from a complex number.
  392.  */
  393. COMPLEX *
  394. csubq(c, q)
  395.     COMPLEX *c;
  396.     NUMBER *q;
  397. {
  398.     COMPLEX *r;
  399.  
  400.     if (qiszero(q))
  401.         return clink(c);
  402.     r = comalloc();
  403.     r->real = qsub(c->real, q);
  404.     r->imag = qlink(c->imag);
  405.     return r;
  406. }
  407.  
  408.  
  409. /*
  410.  * Shift the components of a complex number left by the specified
  411.  * number of bits.  Negative values shift to the right.
  412.  */
  413. COMPLEX *
  414. cshift(c, n)
  415.     COMPLEX *c;
  416.     long n;
  417. {
  418.     COMPLEX *r;
  419.  
  420.     if (ciszero(c) || (n == 0))
  421.         return clink(c);
  422.     r = comalloc();
  423.     r->real = qshift(c->real, n);
  424.     r->imag = qshift(c->imag, n);
  425.     return r;
  426. }
  427.  
  428.  
  429. /*
  430.  * Scale a complex number by a power of two.
  431.  */
  432. COMPLEX *
  433. cscale(c, n)
  434.     COMPLEX *c;
  435.     long n;
  436. {
  437.     COMPLEX *r;
  438.  
  439.     if (ciszero(c) || (n == 0))
  440.         return clink(c);
  441.     r = comalloc();
  442.     r->real = qscale(c->real, n);
  443.     r->imag = qscale(c->imag, n);
  444.     return r;
  445. }
  446.  
  447.  
  448. /*
  449.  * Multiply a complex number by a real number.
  450.  */
  451. COMPLEX *
  452. cmulq(c, q)
  453.     COMPLEX *c;
  454.     NUMBER *q;
  455. {
  456.     COMPLEX *r;
  457.  
  458.     if (qiszero(q))
  459.         return clink(&_czero_);
  460.     if (qisone(q))
  461.         return clink(c);
  462.     if (qisnegone(q))
  463.         return cneg(c);
  464.     r = comalloc();
  465.     r->real = qmul(c->real, q);
  466.     r->imag = qmul(c->imag, q);
  467.     return r;
  468. }
  469.  
  470.  
  471. /*
  472.  * Divide a complex number by a real number.
  473.  */
  474. COMPLEX *
  475. cdivq(c, q)
  476.     COMPLEX *c;
  477.     NUMBER *q;
  478. {
  479.     COMPLEX *r;
  480.  
  481.     if (qiszero(q))
  482.         math_error("Division by zero");
  483.     if (qisone(q))
  484.         return clink(c);
  485.     if (qisnegone(q))
  486.         return cneg(c);
  487.     r = comalloc();
  488.     r->real = qdiv(c->real, q);
  489.     r->imag = qdiv(c->imag, q);
  490.     return r;
  491. }
  492.  
  493.  
  494. /*
  495.  * Take the integer quotient of a complex number by a real number.
  496.  * This is defined to be the result of doing the quotient for each component.
  497.  */
  498. COMPLEX *
  499. cquoq(c, q)
  500.     COMPLEX *c;
  501.     NUMBER *q;
  502. {
  503.     COMPLEX *r;
  504.  
  505.     if (qiszero(q))
  506.         math_error("Division by zero");
  507.     r = comalloc();
  508.     r->real = qquo(c->real, q);
  509.     r->imag = qquo(c->imag, q);
  510.     return r;
  511. }
  512.  
  513.  
  514. /*
  515.  * Take the modulus of a complex number by a real number.
  516.  * This is defined to be the result of doing the modulo for each component.
  517.  */
  518. COMPLEX *
  519. cmodq(c, q)
  520.     COMPLEX *c;
  521.     NUMBER *q;
  522. {
  523.     COMPLEX *r;
  524.  
  525.     if (qiszero(q))
  526.         math_error("Division by zero");
  527.     r = comalloc();
  528.     r->real = qmod(c->real, q);
  529.     r->imag = qmod(c->imag, q);
  530.     return r;
  531. }
  532.  
  533.  
  534. /*
  535.  * Construct a complex number given the real and imaginary components.
  536.  */
  537. COMPLEX *
  538. qqtoc(q1, q2)
  539.     NUMBER *q1, *q2;
  540. {
  541.     COMPLEX *r;
  542.  
  543.     if (qiszero(q1) && qiszero(q2))
  544.         return clink(&_czero_);
  545.     r = comalloc();
  546.     if (!qiszero(q1))
  547.         r->real = qlink(q1);
  548.     if (!qiszero(q2))
  549.         r->imag = qlink(q2);
  550.     return r;
  551. }
  552.  
  553.  
  554. /*
  555.  * Compare two complex numbers for equality, returning FALSE if they are equal,
  556.  * and TRUE if they differ.
  557.  */
  558. BOOL
  559. ccmp(c1, c2)
  560.     COMPLEX *c1, *c2;
  561. {
  562.     BOOL i;
  563.  
  564.     i = qcmp(c1->real, c2->real);
  565.     if (!i)
  566.         i = qcmp(c1->imag, c2->imag);
  567.     return i;
  568. }
  569.  
  570.  
  571. /*
  572.  * Allocate a new complex number.
  573.  */
  574. COMPLEX *
  575. comalloc()
  576. {
  577.     COMPLEX *r;
  578.  
  579.     r = (COMPLEX *) allocitem(&freelist);
  580.     if (r == NULL)
  581.         math_error("Cannot allocate complex number");
  582.     r->links = 1;
  583.     r->real = qlink(&_qzero_);
  584.     r->imag = qlink(&_qzero_);
  585.     return r;
  586. }
  587.  
  588.  
  589. /*
  590.  * Free a complex number.
  591.  */
  592. void
  593. comfree(c)
  594.     COMPLEX *c;
  595. {
  596.     if (--(c->links) > 0)
  597.         return;
  598.     qfree(c->real);
  599.     qfree(c->imag);
  600.     freeitem(&freelist, (FREEITEM *) c);
  601. }
  602.  
  603. /* END CODE */
  604.